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The increasing need for efficient traffic optimization measures is making reliable, fast, 
and robust methods for traffic simulation more and more important. Apart from the de- 
velopment of cellular automata models of traffic flow, this need has stimulated studies of 
suitable numerical algorithms for the solution of macroscopic traffic equations based on 
partial differential equations. 

The numerical integration of partial differential equations is a particularly difficult task, 
and there is no generally applicable method. In contrast to ordinary differential equations, 
the most natural explicit finite difference methods such as replacing differentials by forward, 
backward, or central finite differences often are numerically unstable,00 even in the limit of 
very small discretizations in space and time. For example, a necessary condition for stability 
is that the maximum speed at which information about the solution can be propagated by 
the numerical method must be always greater than the propagation velocity of the exact 
solution.Sii 

In general, the numerical solution of partial differential equations requires special meth- 
ods, which work only under certain conditions. Implicit integration methods are usually 
more stable, but they require the frequent solution of linear systems with multi-diagonal 
matrices. In this column, we will discuss only explicit methods, because they are useful for 
the varying boundary conditions found in realistic traffic simulations, where data is con- 
tinuously fed into the simulation. In addition, explicit methods are more flexible for the 
simulation of on- and off-ramps or entire road networks. 

Macroscopic Traffic Models 

Because the number of vehicles is conserved, all macroscopic traffic models are based on 
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the continuity equation for the vehicle density t) per lane at position x and time t: 

where V{x,t) is the average vehicle velocity. According to Eq. (|I|), the temporal change 
dp/dt of the vehicle density is given by the spatial change —dQ/dx of the traffic flow 
Q{x,t) = p{x,t)V{x,t) and the rate h''^{x,t) > of vehicles entering (+) or leaving (— ) the 
highway at on- or off-ramps. 

To describe time and spatially varying velocities such as occur in emergent traffic jams 
and stop-and-go traffic, we need a dynamical velocity equation. For most continuous models, 
this equation can be written a 

— V— - --— -(V-V) (2) 
dt P '^^ r 

Transport Term Pressure Term Relaxation Term 

According to Eq. (^, the change dV/dt of the average vehicle velocity is given by three 
terms. The transport term originates from the propagation of the velocity profile with the 
velocity V of the vehicles. The pressure term reflects either an anticipation of spatial changes 
in the traffic situation or dispersion effects due to a finite variance of the vehicle velocities. 
The relaxation term describes the adaptation to a dynamic equilibrium velocity Vg with a 
relaxation time r. 

All forms of congested traffic seem to have almost universal properties which are largely 
independent of the initial conditions and the spatially averaged density. For example, the 
characteristic outflow Qout from traffic jams is about 1800 ± 200 vehicles per kilometer and 
lane, and a typical dissolution velocity C of about -15 ± 5 kilometers per hour.i This uni- 
versality arises from the highly correlated state of motion produced by traffic congestions.^ 
In particular, the outflow Qout is related to the time interval between successive departures 
from the traffic jam. Therefore, the outflow is almost independent of the kind and density 
of congested traffic. As a consequence of the constant outflow, the dissolution velocity of 
traffic jams is nearly constant as well. These observations and the transition from free to 



"synchronized" congested trafficEl are correctly described by the nonlocal, gas-kinetic-based 
traffic modehllUll which we now introduce. 



The Nonlocal, Gas-Kinetic-Based Traffic Model 

Our approach is to derive macroscopic traffic models from gas-kinetic traffic 
equations,@0 which are obtained from "microscopic" models of driver-vehicle behavior.@ 
Gas-kinetic traffic models have been proposed earlier, but the correct treatment of the most 
interesting regime of moderate and high densities was still an open problem. 

We have managed to evaluate the vehicular interaction term of the gas-kinetic traffic 
model almost exactly.0 The analytical result can be represented by the nonlocal, dynamical 
equilibrium velocity V^: 



^ d + 0' ( P'T \\ 
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According to Eq. (^, 14 is given by the desired (maximum) velocity Vq, reduced by a term 
which reflects necessary deceleration maneuvers. Here, Pmax is the maximum vehicle density, 
and T is the average time headway at large densities. For example, German authorities 
require that the distance in meters to the front vehicle be not less than half the velocity in km 
per hour, which gives a time headway of T = 1.8 s. For the intra-lane variance 6, we assume 
the constitutive relation 6 = A{p)V^ where A{p) is given by Eq. @. The prime indicates 
that the variable is calculated at the advanced "interaction point" x' = x -l-7(l/pmax + VT) 
with 1 < 7 < 2 rather than at the actual position x. This factor accounts for the fact that 
drivers anticipate the behavior of vehicles in front of them. The monotonically increasing 
"Boltzmann factor" 
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describes^ the dependence of the braking interaction on the dimensionless velocity difference 



Sv = (y — V')/y/6 + 6' between the actual location x and x'. In homogeneous traffic, we 
have -B(O) = 1. If the preceding cars are much slower [Sy ^ 0), the interaction strength 



given by B is particularly large, and it follows that B{Sv) = 'iSy. For Sy 0, we have 
B{Sv) ~ 0. That is, because the distance to the next vehicle is increasing, the vehicle will 
not brake, even if its distance is smaller than the safe distance. 

Finally, the dynamics of the intra-lane variance can be approximated by the constitutive 
relation 



e{x,t) = A{p)V\x,t) 



Aq + AA tanh 



V\x,t), (5) 



Ap 

where the coefficients Aq = 0.008, AA = 0.02, pc = 0.27pniax, and Ap = O.OSpmax have been 
obtained from single-vehicle data on a section of the Dutch motorway A9.0 Simulations 
of sections of other motorways, for example, the German motorway A8, give a somewhat 
lower value of AA = 0.01, which we will use in the following. Equation shows that the 
standard deviation \/6 of vehicle velocities is proportional to the average velocity V, with a 
density-dependent proportionality factor which is small for the density range found in free 
traffic. The velocity variance also enters the gas-kinetic traffic pressure 

P{x,t) = p{x,t)e{x,t) . (6) 

If not explicitly stated otherwise, the simulation results presented here were calculated with 
the model parameters Vq = llOkm/h, r = 32 s, T = 1.8 s, pmax = 160 vehicles/km, and 
7 = 1.2. 

The main difference between the gas-kinetic-based traffic model and other macroscopic 
traffic models is the nonlocal character of the braking term. The nonlocal term in Eq. 
has smoothing properties similar to those of a viscosity term, but its effect is forwardly 
directed and, therefore, more realistic. In contrast, models with an explicit viscosity termB 
lead to unphysical humps in the vehicle density (see Figure p and even negative velocities^ 
(see Figure Our model also has favorable properties with respect to numerical stability 
and integration speed, and hence allows a robust real time simulation of freeway stretches 
up to several thousand kilometers on a personal computer. 

Some Explicit Finite Difference Methods 



A desirable property of the above traffic equations is that they can be formulated in 
terms of a conservation equation with a sink/source term s: 

This form makes it possible to use a variety of numerical standard methods developed for the 
simulation of hydrodynamic problems.0flii The conservative form of the traffic equations 
reads 

and 

dQ ^ c}(gvp + p) _ - g ^ 

dt dx T ' 

where Q{x,t) = p{x,t)V{x,t) is the traffic flow and denotes the average velocity of 
vehicles which enter (+) or leave (— ) the freeway at ramps. For the case ^ V, we obtain 
the additional terms ^u^CV^ — V^)/p in (0). We have 

u=[p,Q], (10) 

f = [Q,{QVp + P)], (11) 

s=[±iy^,{pV,-Q)/r±iy^V^]. (12) 

For the explicit numerical solution methods we will discuss, x and t are discretized with 
uniform values of Aa; and At, respectively. Hence, we calculate u at the discrete points 
(j Aa;, n At) with j, n G {0, 1, 2, . . .}. For brevity, we use the notation = u{j Ax, n At). 
We will discuss the following numerical integration methods:0@ 
1. Lax-Friedrichs method 



^ ^ 2Ax •'^-'^ ^ 
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2. Upwind method 
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3. MacCormack method 
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Consistency Order. If the spatial variation of u is sufficiently smooth, the Lax-Friedrichs 
and upwind methods are first-order, that is, the upper bound of the local error is proportional 
to a if Ax and At are simultaneously decreased by a factor of a. In general, the upwind 
method is not stable and the differential operator needs to be decomposed into parts which 
are treated by upwind and downwind differencing, respectively (the Godunov methodSi^). 
Fortunately, for traffic equations of the form (^, the upwind method is stable (and equivalent 
to the Godunov method). 

The two-step MacCormack and Lax- Wendroff methods are second-order, that is, the 
upper bound of the local error is proportional to a^. However, for shock-like solutions, the 
order is lower for all of the above integration methodsS Note that the predictor step of 
the MacCormack method is an ordinary upwind step, and the corrector step consists of the 
average between the predictor and a "downwind" step with / calculated with the values of 
the predictor. Interchanging the order of the upwind and downwind differencings of the two 
steps has nearly no effect for the equations investigated here. 

Accuracy. Although the discretization errors associated with nonlinear equations are 
difficult to determine, good estimates are usually obtained by doing a local linearization, at 
least for smooth solutions.il Because realistic traffic models can produce sharp gradients, but 



no real shock fronts, the hnearization is apphcable to the numerical treatment of macroscopic 
traffic equations. This analysis shows that the upwind method is more accurate than the 
Lax-Friedrichs method. 

Although the main discretization error of the Lax-Friedrich and upwind methods yields 
a numerical diffusion, which causes a smoothening of shock fronts (see Figures and (b)), 
for the MacCormack and Lax-Wendroff methods there is numerical dispersion, which is asso- 
ciated with a slower propagation of waves with small amplitudes and can lead to oscillations 
in the density and flow fields behind (but not before) large gradients (cf. Figures |^c and 
(d)). However, nonlinear instabilities (see below) also may lead to oscillations. An example 
is the oscillations at the downstream front of the large amplitude jam of Figure 

Numerical Stability. An integration method is numerically unstable if errors grow ex- 
ponentially, which usually leads to wildly oscillating density profiles with very short wave 
lengths and eventually to overfiow error. Even in the quasi-linear case, the above explicit 
discretization methods can lead to three types of instabilities which implies three conditions 
for numerical stability: 

1. Convective instability. Instabilities of the finite difference method for the fiux term 
lead to the Courant-Friedrichs-Lewy condition0ii which for our model becomes 

|At|<^, (19) 

where Vq is the maximum average velocity. For example, for Vq = 40m/s (144km/h) 
and a spatial discretization of Ax = 20 m, we obtain At < 0.5 s. 

2. Diffusion instability. Models which contain an explicit viscosity term Dd^V/dx'^ , 
change the otherwise hyperbolic character of the partial differential equations to a 
parabolic one. For numerical stability, the additional diffusional Courant-Friedrichs- 
Lewy condition D{AxY/{2At) < 1 must be fulfilled.0'ii This condition does not 
apply to simulations of the nonlocal, gas-kinetic-based traffic model. 
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3. Relaxational instability. Due to the finite At, instabilities can develop for spatially 
homogeneous density and flow fields if At is larger than one of the local relaxation 
times 1/rfc, where are the eigenvalues of the functional matrix of the sink/source 
term s. This condition puts further restrictions on At, which depend on the maximum 
vehicle density in the simulation. 

Because viscosity and diffusion terms are replaced by a nonlocal term, coarser discretiza- 
tion is possible in the above gas-kinetic-based traffic model than in most other traffic models, 
allowing real time simulation of much larger systems. For simulations of this model, good 
results are obtained for Ax = 20 m and At = 0.4 s for the MacCormack as well as the upwind 
method. Figure ^ shows that a finer discretization gives almost identical results. 

The development of traffic instabilities starting with almost homogeneous initial traffic is 
a very strict test of numerical accuracy (see Figure In most situations, for example, when 
simulating fronts or already developed congested traffic, the accumulated discretization error 
is much smaller. 

In addition to the quasi-linear instabilities discussed above, genuine nonlinear instabilities 
may arise for certain numerical methods and simulation conditions. An example is the 
oscillations at the downstream front of the large amplitude jam of Figure |^d. It turns out 
that second-order methods are more sensitive to these instabilities. Note that, although 
the quasi-linear behavior is the same for the MacCormack and Lax-Wendroff methods, the 
nonlinear behavior can be different .0 For the above traffic equations, the MacCormack 
method is more stable than the Lax- Wendroff method, and hence we will use the MacCormack 
method whenever a second-order method is desired. 

When simulating the nonlocal model with one of the first-order methods, we did not 
observe nonlinear instabilities. Furthermore, higher order methods are not necessarily more 
accurate^ especially for large gradients. Thus, one should always implement different nu- 
merical methods and compare their simulation results. Aside from their double integration 
speed, it is sometimes preferable to use first-order methods (causing somewhat smoothed 
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wave fronts) instead of second-order methods (producing oscillations and sometimes non- 
linear numerical instabilities). One reason is that oscillations act like perturbations, which 
may give rise to additional traffic jams that do not correspond with reality. 

Initial and Boundary Conditions 

To calculate traffic flows, it is necessary to specify the initial and boundary conditions. 
The initial conditions are completely determined by specifying u{x, 0) in a range [0, L + 5], 
where L is the length of the simulated section and d is the maximum nonlocality, which 
is of the order of 30 m in the gas-kinetic-based model. In systems with open rather than 
periodic boundaries, the initial conditions usually influence the simulation results only for 
a short time, at least, if we starts with free traffic, because errors in the specification of the 
initial conditions are propagated outside the simulated freeway stretch very quickly. Thus, 
the choice of the initial conditions is unimportant. For example, we can start with a linear 
interpolation of the measured initial boundary values. However, because of conservation of 
vehicle number, the initial conditions are relevant for closed systems with periodic boundary 
conditions, which are given by u{0,t) — u{L,t), and du{0,t)/dx — du{L,t)/dx. 

The specification of time-dependent boundary conditions is much more involved. The 
following options are reasonable in different situations: 

1. Dirichlet boundary conditions are given by the empirically measured values w(0, t) and 
u{L, t) at both ends of the particular freeway stretch. For the values at x e [L,L + 5] 
beyond the right boundary (required by the nonlocal term of the gas-kinetic-based 
model), a constant u{x) — u{L) of the boundary values is assumed. 

2. Homogeneous von Neumann boundary conditions assume that the density and vehicular 
flow remain unchanged at the boundaries x — and x — L: 



du{0,t) 



du{L,t) 
dx 



= 0. 



(20) 



dx 



Again, a constant value of u is assumed for x > L. 



9 



3. Free boundary conditions assume that the traffic state is smooth at the boundary, 

9MM)^9MM = o, (21) 

with a linear extrapolation u{L + 5x) = u{L) + 6x du{L)/dx for 6x G [0, 6]. 

4. In- and outflows iQi-mp at on-ramps (+) and off-ramps (— ), can be considered as 
foUows.liHi^ If n is the number of freeway lanes (without ramps) and L^^ap is the 
length of the ramp, we simply set 

i/± = . (22) 

Generally speaking, periodic boundary conditions are best suited for theoretical investi- 
gations of stability, and Dirichlet boundary conditions are best for simulations of real traffic 
with measured values of velocity and the traffic flow at the boundary. Homogeneous von 
Neumann or free boundary conditions for "absorbing" boundaries should be used if the traf- 
fic situation outside the boundaries is not of interest. The effects of the latter two boundary 
conditions are nearly the same. An alternative to using open boundaries is to apply periodic 
boundary conditions for distances which are much longer than the freeway stretch of interest. 
However, for a simulated time interval Tsim = 2h, the required additional length is order 
Vo^^sim ~ 200 km, which considerably reduces the efficiency of the numerical integration. 

There are problems associated with using Dirichlet boundary conditions. Imposing 
Dirichlet boundary conditions at both sides usually leads to an overdetermined system, 
so that either numerical instabilities occur or the boundary conditions are simply ignored 
by the integration method. Even Dirichlet boundary conditions at one of the two bound- 
aries will lead to an unphysical, ill-posed problem in certain situations. This case becomes 
clear when imposing downstream Dirichlet boundary conditions for free traffic flow at low 
densities. If the imposed boundary flow is higher than the flow arriving from the simulated 
section at the boundary, the continuity equation will lead to a decrease of the density, which 
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eventually results in unphysical negative densities. On the other hand, if Dirichlet bound- 
ary conditions are imposed at the upstream boundary in a situation of congested traffic 
where the boundary flow is higher than the equilibrium flow, the continuity equation leads 
to an increased density implying an even lower equilibrium flow. This positive feedback will 
eventually lead to a divergence of the density near the boundary. 

We solve this problem by dynamically switching between Dirichlet and von Neumann 
boundary conditions, depending on the density close to the boundaries. The idea is to 
use Dirichlet boundary conditions if the direction of information flow points towards the 
simulation stretch, and von Neumann boundary conditions otherwise. Because the traffic 
equations contain source terms like the relaxation term, the propagation direction of in- 
formation cannot be obtained from the characteristic velocities (which are the propagation 
velocities of locally linearized homogeneous partial differential equations like conservation 
equations without source terms. Rather, the linearized group velocity Vg — dQ^/ dp of kine- 
matic waves gives a good criterion for separating the traffic situation into free traffic (vg > 0) 
whose information flow points downstream, and congested traffic (vg < 0) where the infor- 
mation flow points upstream. Based on this idea, the following hybrid boundary conditions 
yield good results: If the upstream boundary values po(^) and Qoit) satisfy 

Po{t)<Pip,n or Qo{t) < P2Q{/^x,t) (23) 

Dirichlet boundary conditions are used; pm is the density associated with maximum equi- 
librium flow, /3i = 0.95, and P2 = 0.98. Otherwise, homogeneous Von-Neumann boundary 
conditions are applied. For the downstream boundary conditions, the directions of the in- 
equalities are exchanged. Note that these boundary conditions include situations where 
Dirichlet boundary conditions are necessary at both sides (for example, when a jam enters 
the downstream boundary of the simulation region) as well as situations where no Dirichlet 
boundary conditions are allowed at all (when congested traffic formed at an inhomogeneity 
within the simulated region reaches the upstream boundary). 

We illustrate the effects of different boundary conditions by a simulation of the German 
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motorway A8 near Munich with real traffic data using the upwind method. Figure |^a shows 
a space-time plot of the density with Dirichlet boundary conditions on both sides. The wing- 
like regions of higher density represent "synchronized" congested traffic0 forming upstream 
of a flow-reducing inhomogeneity,^ while the region of the highest density represents a traffic 
jam entering through the downstream boundary. Note that the unnecessary downstream 
boundary conditions for t < 17:30 h are ignored by the upwind differencing, whereas the 
traffic jam entering the downstream boundary for 17:30 h< t < 18:20 h is accepted. This 
desirable behavior cannot be expected in general. For example, the MacCormack method is 
unstable in this situation. 

Figure demonstrates the remarkable property that the boundary conditions for the 
density are nearly irrelevant. Although the boundary conditions for the flow are the same as 
in Fig. a constant amount of 20 vehicles/km has been added artificially to the upstream 
boundary conditions for the density. Nevertheless, the traffic dynamics, in particular the 
spontaneous breakdown at x ~ 40 km and t ~ 17:00 h, is nearly unchanged. This behavior 
can be explained by recognizing that the inflow into the freeway and the outflow from it 
determine the number of vehicles on the freeway stretch, which cannot be changed by the 
traffic dynamics between the boundaries. Hence, changes of the traffic flow are far reaching, 
while the influence of density (or velocity) boundary conditions is restricted to the distance 
that the vehicles drive during the relaxation time. 

Figure |^c illustrates the effect of using homogeneous von Neumann downstream bound- 
ary conditions instead of Dirichlet boundary conditions. For t < 18:00 h, the situation is 
identical to Figure |^a. For t > 18:00 h, however, the upstream moving traffic jam passing 
the downstream boundary is ignored. Hence, the simulation result is very different when 
the direction of the downstream information flow changes. 

Finally, Fig. ^ shows a simulation with the hybrid boundary conditions (see Eq. ([23|)) 
and the upstream boundary shifted downstream by 3.5 km. In this case, the traffic jam 
reaches the upstream boundary, which is handled by the hybrid boundary conditions, while 



12 



all other boundary conditions would lead to a numerical instability for any integration 
method. 

Summary 

We have seen that traffic flows are characterized by the occurrence of congested regions, 
jams, and stop-and-go waves which are associated with large gradients of the vehicle density 
and average velocity. This rich behavior has stimulated an intense research activity. We have 
discussed the problems of macroscopic traffic simulations which are related to the numerical 
integration of systems of coupled nonlinear partial differential equations. Despite the differ- 
ences with hydrodynamic equations, many integration methods developed for conservation 
equations turn out to be apphcable. Compared to implicit methods,0 exphcit methods are 
less robust, but much more flexible with regard to time- dependent boundary conditions and 
(variational) optimization problems, and are usually computationally faster. Among the 
explicit ffist-order methods, the upwind method is more accurate than the Lax-Friedrichs 
method. The MacCormack and the two-step Lax-Wendroff methods are second-order, but 
they are less efficient by a factor of two and produce unrealistic oscillations close to steep 
gradients. The numerical precision may be improved by "high resolution methods," @ where 
a ffist-order method is used for large gradients and a second-order method is used for small 
gradients. This approach may combine the accuracy of the second-order methods with the 
smoothness of the ffist-order methods. 

The most important factor that determines the computation time is the choice of the 
traffic model. In particular, the simulation of the nonlocal, gas-kinetic-based traffic model is 
significantly more efficient than the numerical solution of models with viscosity or diffusion 
terms. This efficiency is mainly related to the fact that the diffusional Courant-Friedrichs- 
Lewy condition, which does not apply for nonlocal terms, is usually far more restrictive than 
the other instability conditions, especially for fine discretizations. Because diffusion terms 
also produce unrealistic effects close to steep gradients, it may be reasonable to generally 
replace models with diffusion or viscosity terms by nonlocal models. Anyway, diffusion and 
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viscosity terms are often a lowest-order approximation of nonlocal terms. They are mainly 
used for historical reasons, because they can be better treated analytically. 

Finally, we have discussed suitable specifications of the boundary conditions. In most 
previous simulation studies, periodic boundary conditions were used to circumvent the in- 
tricate problems related to open boundaries, which are required for the simulation of real 
freeways. We found that Dirichlet boundary conditions work in some cases, but fail in others, 
because of overspecification. The most successful treatment is based on hybrid boundary 
conditions which switch between Dirichlet and homogeneous von Neumann (or free) bound- 
ary conditions depending on the respective direction of information flow. 
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Suggestions for Further Study 

1. Parameter- dependence of flow- density relation, (a) Find the density dependence of the 
stationary, homogeneous solution for the average velocity in the nonlocal, gas-kinetic- 
based traffic model by setting the temporal and spatial derivatives to zero and solving for 
V . (b) The flow-density relation Qe{p) = pVdp) depend on which of the model parameters, 
Vq, r, 7, T, and Pmax? Which parameters are irrelevant for the stationary homogeneous 
solution? (c) Plot Qc and the velocity-density relation Ve for different parameter sets. Start 
with the values used in this column. Check that a speed limit (that is, a decrease in 
Vo) reduces the flow at small densities, but only by a negligible amount for medium and 
high densities. In addition, compare pure car traffic with Vq = 130km/h, T = 1.2 s, and 
Pmax = 160 vehicles/km with traffic containing a considerable fraction of trucks with, for 
example, Vq = 90km/h, T = 3s, and pmax = 110 vehicles/km. In which density regimes do 
the flow-density relations almost agree, and in which are they very different?ffl 
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2. Stability of homogeneous traffic with respect to a localized perturbation. Simulate freeway 
traffic for a circular road of 10 km circumference with the nonlocal, gas-kinetic-based traffic 
model and the parameters used in this column. Assume homogeneous equilibrium traffic 
of density p, and add to the density a localized perturbation of amplitude Ap so that the 
initial conditions are given byii 

pix, 0) = p + Ap [cosh- (^-^] - ^ cosh- 
Q{x,0) = QM- 

Here, xo and xq + Axq are the positions of the positive and negative peaks of the perturbation 
with widths w"^ and w~ , respectively. Choose = 200 m, = 800 m, and Axq = 
w'^ + w~ = 1000 m. Use linear interpolation to calculate the nonlocal terms and take into 
account that the circular road implies that u{L + 6x) = u{6x). Perform the simulations with 
the upwind method on a fixed grid with spatial and temporal discretizations of Ax = 20 m 
and At = 0.4 s. 

(a) Choose a small perturbation amplitude of Ap = 1 vehicle/km and run the simulation 
for 30 minutes. Check that the perturbation does not grow for p < and p > pc3 with 
Pc2 = 29 vehicles/km and pcs = 47 vehicles/km, but gives rise to large-amplitude stop-and- 
go waves for Pc2 < P < Pcs- 

(b) Use larger perturbations of amplitudes up to Ap = 60 vehicles/km and show the 
traffic fiow is metastable. To do so, show that the model produces a single traffic jam 
for p G (pci,Pc2) with pd = 27 vehicles/km and a dipole-hke localized structure for p e 
(pc3,Pc4) with pc4 = 50 vehicles/km, if Ap exceeds a (density-dependent) critical amplitude; 
homogeneous traffic fiow is found for subcritical perturbation amplitudes.0'0 

(c) In sufficiently large systems, there exists a subset of densities p G (pcv,Pc3) in the 
linearly unstable regime where traffic is convectively stable.0 This stability means that the 
localized perturbation will disappear for t — oo at any given location x while, nevertheless, 
the global maximum of the perturbation grows (because the system is linearly unstable). In 
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our case, pcv is given by the value of p at which the downstream boundary of the perturbed 
zone does not move. Illustrate the case of convective stability by making a three-dimensional 
plot of p{x,t) of a simulation running for 15 minutes with p = 45 vehicles/km and Ap = 1 
vehicles/km. Show that pcv = 42.5 vehicles/km. 

3. Characteristic parameters^ Assume a circular road, model parameters, and initial con- 
ditions as in Problem 2. 

(a) Evaluate the densities of fully developed traffic jams and of free traffic at the outflows 
of jams. Determine the associated flows and the propagation velocity of jam fronts. Hint: 
To obtain the densities and flows, run the simulation for 60 minutes and analyze the flelds 
p{x,t) and Q{x,t). The jam density (density of free traffic) is simply the global maximum 
(minimum) of the density, and the jam flow (outflow from jams) is the global minimum 
(maximum) of the flows. To obtain the group velocity, save the density as a function of 
time at a given position (for example, x = 10 km) and determine the time intervals needed 
for the jams to propagate around the circular road. Use the last interval to avoid effects of 
transients. 

(b) Plot the densities, the flows, and the propagation velocity in separate plots as a 
function of the average density p for a perturbation amplitude Ap = 10 vehicle/km. What 
do you flnd? Check that there is a density range in which the density inside and outside of 
traffic jams, the associated flows, and the propagation velocity are independent of p, so that 
these quantities are characteristic parameters of traffic flows. 

(c) Show numerically that the traffic inside and outside of jams is nearly in equilibrium. 
Also show that the group velocity can be expressed analytically in terms of the densities 
and flows inside and outside of jams. This relation means that only two of the characteristic 
parameters, for example, the group velocity and the outflow from jams, are independent 
quantities. 

(d) Do Problems 2a, 2b, and 3a for different values of the relaxation time r. Show 
that with increasing r, the region of linearly unstable traffic and the regions of metastable 
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traffic and the jam density increase, while the outflow from jams decreases. In particular, 
homogeneous traffic is always stable for r < 18 s. 

4. ■■Synchrom.er congested traffiM This trrfBc state occurs it a perturbation reaches a 

stationary inhomogeneity of the road which can be an on-ramp for example. 

(a) Simulate a freeway of length 10 km with open boundaries and an on-ramp of length 
Lrrap = 400 m at Xrmp = 5 km by assuming in Eqs. (||) and (||) a source term z/+ = Qrmp/ L 
inside the ramp region x G [xrmp — -^/2,Xrmp + L/2], but = otherwise. Assume 
r = 40 s, the usual values for the other model parameters, and homogeneous equilibrium 
traffic of density p = 15 vehicles/km without perturbations as initial conditions. To obtain 
stationary conditions, simulate the ffist 20 minutes with a constant ramp flow Qrmp = 
500 vehicles/h using homogeneous von Neumann boundary conditions at both sides of the 
road. Now, introduce a perturbation of amplitude AQrmp = 150 vehicles/km by linearly 
increasing the ramp flow up to Qrmp + ^Qrmp at t = 22.5 min, and decreasing it again to 

min. After running the simulation for a total of 60 minutes, you should 
see "synchronized" congested traffic, that is, an increasing region of high density and low 
velocity, but relatively high flow, whose upstream front is propagating upstream, while the 
downstream front is pinned at the ramp. 

(b) Verify that "synchronized" congested traffic is in equilibrium and determine the 
numerical value of its outflow Qout- Show that the propagation velocity Vg of the upstream 
front can be expressed by the relation 

Qout Qrmp Qmain 
Vg = , 

Pcong(Qout Qrmp) Pfrec(Qmain) 

where Qrmp and Qmain = Qe(p) are the inflows at the ramp and to the main road. The 
densities Pcong(Q) and pfree(Q) denote, in accordance with Qe = pK(p), congested traffic 
(p > 31 vehicles/km) and free traffic (p < 31 vehicles/km), respectively. 

5. Different traffic states close to an on-ramp^ Varying the initial traffic density, the 
length of the on-ramp, and the ramp flow leads to several interesting states of congested 
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traffic. Do simulations as in Problem 4 with (a) p — 20 vehicles/km, Qrmp = 250 vehicles/h, 
AQrmp = 150 vehicles/km, (b) p = 20 vehicles/km, Qrmp = 150 vehicles/h, AQrmp = 150 
vehicles/km, and (c) p = 17 vehicles/km, Qrmp = 250 vehicles/h, AQ^ap = 550 vehicles/km. 
You should observe (a) oscillating congested traffic, (b) triggered stop-and-go waves, and (c) 
pinned localized congestion. By varying the inflow Qmain = Qe{p) to the main road, verify 
that oscillating congested traffic emerges, if the expression for Vg in Problem 4 is negative 
(upstream moving front); otherwise pinned localized congestion or free traffic occur. 
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FIG. 1. Simulation of discontinuous initial conditions for (a) the non-local, gas-kinetic-based 
traffic model, and (b) for a typical traffic model with a viscosity termS The density profile is 

shown at t = Os (— ), 10s (- -), 30s (- - -), 60s (• • •), 120s ( ), and 240s (- • •)• The inset of 

(b) shows an unrealistic detail of the profile (upstream front) at t = 10s. Also note that, in (b), 
the upstream front propagates too fast. 
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FIG. 2. First stages of the density and velocity profiles evolving from a discontinuous upstream 
front (solid lines) corresponding, for example, to an end of a traffic jam behind a curve. The 

lines correspond to t = Os (— ), 5s (- -), 10s ( ), 20s (•••), 60s (- • ~), and 120s (- • •)• 

Parts (a) and (b) show simulations with the non-local, gas-kinetic-based model (Ax = 20 m, 
At = 0.1 s), while (c) and (d) correspond to a typical traffic model with a viscosity termJH§ To obtain 
comparable equilibrium velocities in the initial jam, the initial jam densities /Oj^^^ = 140 vehicles/h 



Ac/d) 



and pj^'^' = 105 vehicles/h were chosen differently. In the viscosity traffic model , the finite velocity 
diffusion leads to an increase of velocity also in the congested part (see (d)) and to a subsequent 
further increase of density (see (c)). The ensuing positive density gradient would lead to negative 
velocities for an initial jam density higher than 106 vehicles/km, which is much lower than the 
maximum density pmax = 160 vehicles/km. 
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FIG. 3. Comparison of the upwind method, (a) and (c), with the MacCormack method, 
(b) and (d), for simulations of the nonlocal, gas-kinetic-based traffic model with discontinuous 
initial conditions ( — ). The plots correspond to transitions between free traffic at two different 
vehicle densities (see (a) and (b)), and to transitions to and from congested traffic with backwards 

propagating fronts (see (c) and (d)). The density is shown at t = Os ( — ), 10s ( — ), 30s ( ), 60s 

(• • •), and 120 s (- • -) (only in (c) and (d)). In the MacCormack method simulations, the oscillations 
behind the large gradient result from the dispersion error (see (b)), whereas the oscillations around 
X = 3 km originate from nonlinear instabilities (see (d)). 
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FIG. 4. Spatio-temporal formation of stop-and-go waves, simulated with the non-local, 
gas-kinetic-based traffic model, using different integration methods and discretizations, (a) Mac- 
Cormack method with At = 0.4 s and Ax = 20 m. (b) Upwind method with At = 0.4 s and 
Ax = 20 m. (c) Upwind method with At = 0.01s and Ax = 20 m. (d) Upwind method with 
At = 0.1 s and Ax = 5 m. 
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FIG. 5. Simulation with empirical boundary conditions at the German freeway A8 near Mu- 
nich, with the model parameters Vq = 130 km/h for x < 41 km, Vq = 97 km/h for x > 41 km, 
r = 32 s, T = 2.0 s, pmax = HO vehicles/km, and 7 = 1.2, ^0 = 0.008, = 0.01, pc = 0.27pmax, 
and Ap = O.lpmax- Note that the decreased desired velocity for x > 41km reflects a flow-reducing 
gradient on the freeway stretch, (a) Dirichlet boundary conditions, applied to p and Q at the up- 
stream and downstream boundaries. We can distinguish two different kinds of congestion: a region 
of "synchronized" trafficS'Ei with increased homogeneous density, and a backwards propagating 
traffic jam entering the freeway stretch at the downstream boundary. Dirichlet boundary condi- 
tions are applicable here, because the traffic jam does not reach the upstream boundary. (The jam 
dissolves before, as the inflow decreases in the course of time.) (b) As in (a), but with a constant 
Ap = 20 vehicles/km added to the measured density at the upstream boundary. The congested 
traffic enforced at the boundary relaxes to free traffic very quickly. Note that the downstream 
traffic patterns do not change significantly, because the boundary flows (and, hence, the number 
of entering and leaving vehicles per time unit) are the same as in (a). Nevertheless, this result 
suggests that the congested traffic patterns are self-organized structures, (c) As in (a), but with 
homogeneous von Neumann boundary conditions for flow and density at the downstream boundary, 
which ignore the upstream moving traffic jam entering at the downstream boundary. However, the 
homogeneous congested region of "synchronized" traffic is correctly reproduced, which indicates 
that it is not triggered by the downstream boundary condition, but it is rather generated by the 
inflow to the motorway at the upstream boundary together with the flow-reducing gradient begin- 
ning at X = 41.0km. (d) Simulation with the hybrid boundary conditions (p3[). Compared to (a), 
we shifted the upstream boundary by 3.5 km in downstream direction, so that it is reached by the 
jam. The resulting spatio-temporal traffic patterns are nearly identical as in (a), indicating that 
the hybrid boundary conditions handle dynamically emerging congestion in a natural way, even at 
the boundaries. 
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